Celltype Colours¶

In [1]:
seurat_clusters_colors <- c(
  "0"  = "#FF6666",   # Red
  "1"  = "#FF9933",   # Orange
  "2"  = "#FFCC99",   # Light Peach
  "3"  = "#FFFF99",   # Light Yellow
  "4"  = "#CCFF99",   # Light Green
  "5"  = "#99FFCC",   # Light Aqua
  "6"  = "#66CCCC",   # Light Cyan
  "7"  = "#66CC99",   # Light Mint
  "8"  = "#66CC66",   # Light Green
  "9"  = "#66CCFF",   # Light Blue
  "10" = "#6666FF",   # Blue
  "11" = "#9966FF",   # Purple
  "12" = "#CC66FF",   # Lavender
  "13" = "#FF6699",   # Pink
  "14" = "#FF9966",   # Orange
  "15" = "#FFCC66",   # Yellow
  "16" = "#FFFF66",   # Yellow
  "17" = "#FF6666",   # Red
  "18" = "#6600FF",   # Indigo
  "19" = "#00FF00",   # Green
  "20" = "#FF00FF",   # Magenta
  "21" = "#00FFFF",   # Cyan
  "22" = "#0000FF",   # Blue
  "23" = "#FFFF00",   # Yellow
  "24" = "#FF0000",   # Red
  "25" = "#800080",   # Purple
  "26" = "#800000",   # Maroon
  "27" = "#808000",   # Olive
  "28" = "#008080",   # Teal
  "29" = "#808080",   # Gray
  "30" = "#C0C0C0",   # Silver
  "31" = "#FFD700",   # Gold
  "32" = "#FF4500",   # OrangeRed
  "33" = "#DA70D6",   # Orchid
  "34" = "#EEE8AA",   # PaleGoldenRod
  "35" = "#98FB98",   # PaleGreen
  "36" = "#AFEEEE",   # PaleTurquoise
  "37" = "#DB7093",   # PaleVioletRed
  "38" = "#FFEFD5",   # PapayaWhip
  "39" = "#FFDAB9",   # PeachPuff
  "40" = "#CD853F",   # Peru
  "41" = "#FFC0CB",   # Pink
  "42" = "#DDA0DD",   # Plum
  "43" = "#B0E0E6",   # PowderBlue
  "44" = "#BC8F8F",   # RosyBrown
  "45" = "#4169E1",   # RoyalBlue
  "46" = "#8B4513",   # SaddleBrown
  "47" = "#FA8072",   # Salmon
  "48" = "#FAA460",   # SandyBrown
  "49" = "#2E8B57",   # SeaGreen
  "50" = "#FFF5EE",   # SeaShell
  "51" = "#A0522D",   # Sienna
  "52" = "#C0C0C0",   # Silver
  "53" = "#87CEEB",   # SkyBlue
  "54" = "#6A5ACD"    # SlateBlue
)
In [2]:
subcluster_colors <- c(
  "Epiblast" = "#FF0000",                    # Red
  "Primitive Streak" = "#FF7F00",            # Orange
  "Primitive Streak & PGC" = "#FFFF00",      # Yellow
  "Unresolved Primitive Streak" = "#7FFF00", # Chartreuse
  "Anterior Primitive Streak" = "#00FF00",   # Green
  "Unresolved Mesoderm" = "#00FF7F",         # Spring Green
  "Nascent Mesoderm" = "#00FFFF",            # Cyan
  "Early Haematoendothelial Progenitors" = "#007FFF", # Azure
  "Late Haematoendothelial Progenitors" = "#0000FF",  # Blue
  "Mutant Nascent Mesoderm" = "#7F00FF",     # Violet
  "Unresolved Mutant Mesoderm" = "#FF00FF",  # Magenta
  "Mutant Lateral Plate Mesoderm #1" = "#FF007F", # Rose
  "Mutant Lateral Plate Mesoderm #2" = "#7F0000"  # Dark Red
)
In [3]:
stages_colors = c(
"E6.5" = "#F21A00",
"E6.75" = "#EC4B00",
"E7.0" = "#E67D00",
"E7.25" = "#E1AF00",
"E7.5" = "#E4B80E",
"E7.75"="#E7C21C",
"E8.0" = "#EBCC2A",
"E8.25" = "#C4C55D",
"E8.5" = "#9EBE91",
"E8.75" = "#78B7C5",
"E9.0" = "#3399FF",
"E9.25"="#297ACC",
"E9.5"="#2162A3",
"Mixed gastrulation" = "#BEBEBE")
In [4]:
stages_colors_updated = c(
  "E6.5" = "#FF0000",       # Red
  "E6.75" = "#FF4500",      # Orange-Red
  "E7.0" = "#FFA500",       # Orange
  "E7.25" = "#FFD700",      # Gold
  "E7.5" = "#FFFF00",       # Yellow
  "E7.75" = "#00BFFF",      # Deep Sky Blue
  "E8.0" = "#1E90FF",       # Dodger Blue
  "E8.25" = "#4169E1",      # Royal Blue
  "E8.5" = "#800080"        # Purple
)
In [5]:
anatomy_colors = c(
"EP" = "#000000",
"YS" = "#5E4FA2",
"Anterior section" = "#774611",
"Posterior section" = "#A2D630",
"Medial section" = "#D95F02",
"Posterior"="#E2C207",
"Anterior" = "#F71616",
"Pooled" = "#BEBEBE",
"NA" = "#BEBEBE")
In [6]:
# We can add colours as well if you really want!
# this is the colour scheme, we usually keep this in a seperate file that we then call in our notebooks
celltype.colors = c(
  "Epiblast" = "#635547",
  "Primitive_Streak" = "#DABE99",
  "Caudal_epiblast" = "#9e6762",
  "PGC" = "#FACB12",
  "Anterior_Primitive_Streak" = "#c19f70",
  "Notochord" = "#0F4A9C",
  "Def._endoderm" = "#F397C0",
  "Gut" = "#EF5A9D",
  "Nascent_mesoderm" = "#C594BF",
  "Mixed_mesoderm" = "#DFCDE4",
  "Intermediate_mesoderm" = "#139992",
  "Caudal_Mesoderm" = "#3F84AA",
  "Paraxial_mesoderm" = "#8DB5CE",
  "Somitic_mesoderm" = "#005579",
  "Pharyngeal_mesoderm" = "#C9EBFB",
  "Cardiomyocytes" = "#B51D8D",
  "Allantois" = "#532C8A",
  "ExE_mesoderm" = "#8870ad",
  "Mesenchyme" = "#cc7818",
  "Haematoendothelial_progenitors" = "#FBBE92",
  "Endothelium" = "#ff891c",
  "Blood_progenitors" = "#c9a997",
  "Blood_progenitors_1" = "#f9decf",
  "Blood_progenitors_2" = "#c9a997",
  "Erythroid" = "#EF4E22",
  "Erythroid1" = "#C72228",
  "Erythroid2" = "#f79083",
  "Erythroid3" = "#EF4E22",
  "NMP" = "#8EC792",
  "Neurectoderm" = "#65A83E",
  "Rostral_neurectoderm" = "#65A83E",
  "Caudal_neurectoderm" = "#354E23",
  "Neural_crest" = "#C3C388",
  "Forebrain_Midbrain_Hindbrain" = "#647a4f",
  "Spinal_cord" = "#CDE088",
  "Surface_ectoderm" = "#f7f79e",
  "Visceral_endoderm" = "#F6BFCB",
  "ExE_endoderm" = "#7F6874",
  "ExE_ectoderm" = "#989898",
  "Parietal_endoderm" = "#1A1A1A",
  "Lateral_plate_mesoderm" = "#F9DFE6"
)
In [7]:
original.colors = c(
  "Epiblast" = "#635547",
  "Primitive Streak" = "#DABE99",
  "Anterior Primitive Streak" = "#c19f70",
  "Caudal epiblast" = "#9e6762",
  "PGC" = "#FACB12",
  "Anterior Primitive_Streak" = "#c19f70",
  "Notochord" = "#0F4A9C",
  "Def. endoderm" = "#F397C0",
  "Gut" = "#EF5A9D",
  "Nascent mesoderm" = "#C594BF",
  "Mixed mesoderm" = "#DFCDE4",
  "Intermediate mesoderm" = "#139992",
  "Caudal Mesoderm" = "#3F84AA",
  "Paraxial mesoderm" = "#8DB5CE",
  "Somitic mesoderm" = "#005579",
  "Pharyngeal mesoderm" = "#C9EBFB",
  "Cardiomyocytes" = "#B51D8D",
  "Allantois" = "#532C8A",
  "ExE mesoderm" = "#8870ad",
  "Mesenchyme" = "#cc7818",
  "Haematoendothelial progenitors" = "#FBBE92",
  "Endothelium" = "#ff891c",
  "Blood progenitors" = "#c9a997",
  "Blood progenitors 1" = "#f9decf",
  "Blood progenitors 2" = "#c9a997",
  "Erythroid" = "#EF4E22",
  "Erythroid1" = "#C72228",
  "Erythroid2" = "#f79083",
  "Erythroid3" = "#EF4E22",
  "NMP" = "#8EC792",
  "Neurectoderm" = "#65A83E",
  "Rostral neurectoderm" = "#65A83E",
  "Caudal neurectoderm" = "#354E23",
  "Neural crest" = "#C3C388",
  "Forebrain Midbrain_Hindbrain" = "#647a4f",
  "Spinal cord" = "#CDE088",
  "Surface ectoderm" = "#f7f79e",
  "Visceral endoderm" = "#F6BFCB",
  "ExE endoderm" = "#7F6874",
  "ExE ectoderm" = "#989898",
  "Parietal endoderm" = "#1A1A1A"
)
In [8]:
extended.colors = c(
"Epiblast" = "#635547",
"Primitive Streak" = "#DABE99",
"Caudal epiblast" = "#9e6762",
"PGC" = "#FACB12",
"Anterior Primitive Streak" = "#c19f70",
"Node"="#153b3d",
"Notochord" = "#0F4A9C",
"Gut tube" = "#EF5A9D",
"Hindgut" = "#F397C0",
"Midgut" = "#ff00b2",
"Foregut" = "#ffb7ff",
"Pharyngeal endoderm"="#95e1ff",
"Thyroid primordium"="#97bad3",
"Nascent mesoderm" = "#C594BF",
"Intermediate mesoderm" = "#139992",
"Caudal mesoderm" = "#3F84AA",
"Lateral plate mesoderm" = "#F9DFE6",
"Limb mesoderm" = "#e35f82",
"Forelimb" = "#d02d75",
"Kidney primordium" = "#e85639",
"Presomitic mesoderm"="#5581ca",#"#0000ff",#blue
"Somitic mesoderm" = "#005579",
"Posterior somitic tissues" = "#5adbe4",#"#40e0d0",#turquoise
"Paraxial mesoderm" = "#8DB5CE",
"Cranial mesoderm" = "#456722",#"#006400",#darkgreen
"Anterior somitic tissues"= "#d5e839",
"Sclerotome" = "#e3cb3a",#"#ffff00",#yellow
"Dermomyotome" = "#00BFC4",#"#a52a2a",#brown
"Pharyngeal mesoderm" = "#C9EBFB",
"Cardiopharyngeal progenitors" = "#556789",
"Anterior cardiopharyngeal progenitors"="#683ed8",
"Allantois" = "#532C8A",
"Mesenchyme" = "#cc7818",
"YS mesothelium" = "#ff7f9c",
"Epicardium"="#f79083",
"Embryo proper mesothelium" = "#ff487d",
"Cardiopharyngeal progenitors FHF"="#d780b0",
"Cardiomyocytes FHF 1"="#a64d7e",
"Cardiomyocytes FHF 2"="#B51D8D",
"Cardiopharyngeal progenitors SHF"="#4b7193",
"Cardiomyocytes SHF 1"="#5d70dc",
"Cardiomyocytes SHF 2"="#332c6c",
"Haematoendothelial progenitors" = "#FBBE92",
"Blood progenitors" = "#6c4b4c",
"Erythroid" = "#C72228",
"Chorioallantoic-derived erythroid progenitors"="#E50000",
"Megakaryocyte progenitors"="#e3cb3a",
"MEP"="#EF4E22",
"EMP"="#7c2a47",
"YS endothelium"="#ff891c",
"YS mesothelium-derived endothelial progenitors"="#AE3F3F",
"Allantois endothelium"="#2f4a60",
"Embryo proper endothelium"="#90e3bf",
"Venous endothelium"="#bd3400",
"Endocardium"="#9d0049",
"NMPs/Mesoderm-biased" = "#89c1f5",
"NMPs" = "#8EC792",
"Ectoderm" = "#ff675c",
"Optic vesicle" = "#bd7300",
"Ventral forebrain progenitors"="#a0b689",
"Early dorsal forebrain progenitors"="#0f8073",
"Late dorsal forebrain progenitors"="#7a9941",
"Midbrain/Hindbrain boundary"="#8ab3b5",
"Midbrain progenitors"="#9bf981",
"Dorsal midbrain neurons"="#12ed4c",
"Ventral hindbrain progenitors"="#7e907a",
"Dorsal hindbrain progenitors"="#2c6521",
"Hindbrain floor plate"="#bf9da8",
"Hindbrain neural progenitors"="#59b545",
"Neural tube"="#233629",
"Migratory neural crest"="#4a6798",
"Branchial arch neural crest"="#bd84b0",
"Frontonasal mesenchyme"="#d3b1b1",
"Spinal cord progenitors"="#6b2035",
"Dorsal spinal cord progenitors"="#e273d6",
"Non-neural ectoderm" = "#f7f79e",
"Surface ectoderm" = "#fcff00",
"Epidermis" = "#fff335",
"Limb ectoderm" = "#ffd731",
"Amniotic ectoderm" = "#dbb400",
"Placodal ectoderm" = "#ff5c00",
"Otic placode"="#f1a262",
"Otic neural progenitors"="#00b000",
"Visceral endoderm" = "#F6BFCB",
"ExE endoderm" = "#7F6874",
"ExE ectoderm" = "#989898",
"Parietal endoderm" = "#1A1A1A"
)
In [9]:
integrated_colors_updated <- c(
  "ExE Ectoderm" = "#989898",                      # Red
  "ExE Endoderm" = "#7F6874",                      # Orange
  "Visceral Endoderm" = "#F6BFCB",                 # Yellow
  "Epiblast" = "#635547",                          # Green
  "Primitive Streak" = "#DABE99",                  # Cyan
  "Anterior Primitive Streak" = "#FFA500",         # Orange
  "Nascent mesoderm #1" = "#C594BF",               # Magenta
  "Nascent mesoderm #2" = "#8C5A82",               # Darker Purple
  "Haematoendothelial progenitors #1" = "#FBBE92", # Rose
  "Haematoendothelial progenitors #2" = "#6c4b4c", # Chartreuse
  "Posterior Primitive Streak" = "#5A5AC5",        # Soft Blue
  "Posterior Lateral Plate Mesoderm" = "#7F0000",  # Maroon
  "Mesenchyme" = "#cc7818",                        # Azure
  "Allantois Mesoderm" = "#3D99C3",                # Medium Blue
  "Intermediate Mesoderm" = "#139992",             # Spring Green
  "Low Quality" = "#808080",                       # Gray
  "Parietal Endoderm" = "#1A1A1A"                  # Black
)

Load in Packages¶

In [19]:
suppressPackageStartupMessages({ 
    library(data.table) 
    library(dplyr) 
    library(ggplot2) 
    library(SingleCellExperiment)
    library(dplyr)
    library(celldex)
    library(SingleR)
    library(RColorBrewer)
    library(scater) 
    library(StabMap) 
    library(scran) 
    library(harmony) 
    library(patchwork)
    library(Seurat)
    library(plotly)
    library(Matrix)
    library(stringr)
    library(ggalluvial)
    library(NMF)
    library(CellChat)
    library(patchwork)
    library(ComplexHeatmap)
})

Set Data Locations and Load in the Datasets¶

In [70]:
# set paths to data locations
io = list()
io$main = "/rds/project/rds-SDzz0CATGms/users/ltgh2/"
# set the working directory
setwd(io$main)
In [71]:
# set paths to data locations for the query samples
io$query = file.path(io$main, "projects/10_Axin_1_2/github_submission/outputs/1/integrated_axin_dataset_seurat.rds")

# load in the query data 
query <- readRDS(io$query)
In [72]:
# Load in the new MetaData from the Mapping
# set paths to data locations for the query samples
io$axin_transferred_anatomy = file.path(io$main, "projects/10_Axin_1_2/outputs/axin_transferred_anatomy.rds")
axin_transferred_anatomy <- readRDS(io$axin_transferred_anatomy)

io$axin_transferred_celltype_extended_atlas = file.path(io$main, "projects/10_Axin_1_2/outputs/axin_transferred_celltype_extended_atlas.rds")
axin_transferred_celltype_extended_atlas <- readRDS(io$axin_transferred_celltype_extended_atlas)

io$axin_transferred_celltype_PijuanSala2019 = file.path(io$main, "projects/10_Axin_1_2/outputs/axin_transferred_celltype_PijuanSala2019.rds")
axin_transferred_celltype_PijuanSala2019 <- readRDS(io$axin_transferred_celltype_PijuanSala2019)

io$axin_transferred_stage = file.path(io$main, "projects/10_Axin_1_2/outputs/axin_transferred_stage.rds")
axin_transferred_stage <- readRDS(io$axin_transferred_stage)

io$integrated_object_clusters = file.path(io$main, "projects/10_Axin_1_2/outputs/integrated_object_clusters.rds")
integrated_object_clusters <- readRDS(io$integrated_object_clusters)
In [73]:
# subset the predicted labels and the mapping score
axin_extended_celltype <- subset(axin_transferred_celltype_extended_atlas, select = c(predicted.id, prediction.score.max))
colnames(axin_extended_celltype) <- c("extended_celltype", "extended_celltype_mappingscore")          

axin_original_celltype <- subset(axin_transferred_celltype_PijuanSala2019, select = c(predicted.id, prediction.score.max))
colnames(axin_original_celltype) <- c("original_celltype", "original_celltype_mappingscore")    
                                 
axin_stage <- subset(axin_transferred_stage, select = c(predicted.id, prediction.score.max))
colnames(axin_stage) <- c("stage", "stage_mappingscore")    
                     
axin_anatomy <- subset(axin_transferred_anatomy, select = c(predicted.id, prediction.score.max))
colnames(axin_anatomy) <- c("anatomy", "anatomy_mappingscore")    

# make a complete metadata table with all the transferred labels
axin_transferred_labels <- cbind(axin_extended_celltype, axin_original_celltype, axin_stage, axin_anatomy)
In [74]:
tail(integrated_object_clusters)
A data.frame: 6 x 1
integrated_clusters
<fct>
LibA_DKOAxin1Axin2_E6_1_IGO_14590_1#TTTGGTTAGACGTCGA-129
LibA_DKOAxin1Axin2_E6_1_IGO_14590_1#TTTGGTTGTCTGTCCT-17
LibA_DKOAxin1Axin2_E6_1_IGO_14590_1#TTTGGTTTCCGCGGAT-11
LibA_DKOAxin1Axin2_E6_1_IGO_14590_1#TTTGGTTTCTCAACGA-123
LibA_DKOAxin1Axin2_E6_1_IGO_14590_1#TTTGTTGAGCGCTGCT-122
LibA_DKOAxin1Axin2_E6_1_IGO_14590_1#TTTGTTGGTCTGTGGC-124

Add the MetaData to the Axin Dataset¶

In [75]:
# add the new cell type labels to the axin integrated object
query <- AddMetaData(query, metadata = axin_transferred_labels)
In [76]:
query <- AddMetaData(query, metadata = integrated_object_clusters)
In [77]:
colnames(query@meta.data)
  1. 'orig.ident'
  2. 'nCount_RNA'
  3. 'nFeature_RNA'
  4. 'barcode'
  5. 'sample'
  6. 'cell'
  7. 'mitochondrial_percent_RNA'
  8. 'ribosomal_percent_RNA'
  9. 'idx'
  10. 'stage'
  11. 'tdTom'
  12. 'tdTom_corr'
  13. 'pass_rnaQC'
  14. 'doublet_score'
  15. 'doublet_call'
  16. 'celltype.mapped_mnn'
  17. 'celltype.score_mnn'
  18. 'celltype_extended.mapped_mnn'
  19. 'celltype_extended.score_mnn'
  20. 'stage.mapped_mnn'
  21. 'cellstage.score_mnn'
  22. 'closest.cell_mnn'
  23. 'genotype'
  24. 'RNA_snn_res.1'
  25. 'seurat_clusters'
  26. 'across_genotype_integration_snn_res.1'
  27. 'extended_celltype'
  28. 'extended_celltype_mappingscore'
  29. 'original_celltype'
  30. 'original_celltype_mappingscore'
  31. 'stage_mappingscore'
  32. 'anatomy'
  33. 'anatomy_mappingscore'
  34. 'integrated_clusters'
In [78]:
options(repr.plot.width=15, repr.plot.height=8) 
object <- query

p1 <- DimPlot(object, reduction = "umap",
        #label = TRUE,
        group.by = "extended_celltype"
       )  + theme(legend.position='none') + scale_color_manual(values=extended.colors) + coord_fixed()

p2 <- FeaturePlot(object, reduction = "umap",
        #label = FALSE,
        feature = "extended_celltype_mappingscore",
        cols = c("red", "blue")
       )  + coord_fixed()

p1 | p2

Rename the Integrated Subcluster¶

In [141]:
object@meta.data$integrated_clusters_character <- as.character(object@meta.data$integrated_clusters)
In [142]:
object$integrated_clusters_updated <- recode(object$integrated_clusters_character, 
  "0" = "Epiblast",
  "1" = "ExE Ectoderm",
  "2" = "Epiblast",
  "3" = "Low Quality",
  "4" = "Nascent mesoderm #1",
  "5" = "Low Quality",
  "6" = "ExE Endoderm",
  "7" = "ExE Ectoderm",
  "8" = "Primitive Streak",
  "9" = "Allantois Mesoderm",
  "10" = "Intermediate Mesoderm",
  "11" = "Primitive Streak",
  "12" = "Low Quality",
  "13" = "ExE Endoderm",
  "14" = "Low Quality",
  "15" = "Primitive Streak",
  "16" = "Low Quality",
  "17" = "Low Quality",
  "18" = "Low Quality",
  "19" = "Posterior Lateral Plate Mesoderm",
  "20" = "Anterior Primitive Streak",
  "21" = "Mesenchyme",
  "22" = "Nascent mesoderm #2",
  "23" = "Visceral Endoderm",
  "24" = "Visceral Endoderm",
  "25" = "Low Quality",
  "26" = "Low Quality",
  "27" = "Epiblast",
  "28" = "Low Quality",
  "29" = "Haematoendothelial progenitors #2",
  "30" = "Haematoendothelial progenitors #1",
  "31" = "ExE Ectoderm",
  "32" = "Low Quality",
  "33" = "Low Quality",
  "34" = "Low Quality",
  "35" = "Posterior Primitive Streak",
  "36" = "Low Quality",
  "37" = "Low Quality",
  "38" = "Low Quality",
  "39" = "Low Quality",
  "40" = "Parietal Endoderm",
  "41" = "Low Quality",
  "42" = "Low Quality"
)
In [81]:
table(object@meta.data$integrated_clusters_updated)
               Allantois Mesoderm         Anterior Primitive Streak 
                              523                               158 
                         Epiblast                      ExE Ectoderm 
                             1596                              1374 
                     ExE Endoderm Haematoendothelial progenitors #1 
                             2151                                88 
Haematoendothelial progenitors #2             Intermediate Mesoderm 
                              169                                71 
                      Low Quality                        Mesenchyme 
                               87                                48 
              Nascent mesoderm #1               Nascent mesoderm #2 
                              724                               493 
                Parietal Endoderm  Posterior Lateral Plate Mesoderm 
                               65                               285 
       Posterior Primitive Streak                  Primitive Streak 
                              175                               807 
                Visceral Endoderm 
                             1436 
In [144]:
# Assuming 'object' is your Seurat object containing the updated cluster information
# Define the desired order of levels
cluster_order <- c("ExE Ectoderm", "ExE Endoderm", "Visceral Endoderm", "Parietal Endoderm",
                   "Epiblast", "Primitive Streak", "Anterior Primitive Streak", "Nascent mesoderm #1", "Nascent mesoderm #2",
                   "Haematoendothelial progenitors #1", "Haematoendothelial progenitors #2", 
                   "Posterior Lateral Plate Mesoderm", "Mesenchyme", "Posterior Primitive Streak", "Allantois Mesoderm", "Intermediate Mesoderm", "Low Quality")
In [145]:
# Convert 'seurat_clusters_updated' to a factor with the specified levels and order
object$integrated_clusters_updated_ordered <- factor(object$integrated_clusters_updated, levels = cluster_order)
In [84]:
colnames(object@meta.data)
  1. 'orig.ident'
  2. 'nCount_RNA'
  3. 'nFeature_RNA'
  4. 'barcode'
  5. 'sample'
  6. 'cell'
  7. 'mitochondrial_percent_RNA'
  8. 'ribosomal_percent_RNA'
  9. 'idx'
  10. 'stage'
  11. 'tdTom'
  12. 'tdTom_corr'
  13. 'pass_rnaQC'
  14. 'doublet_score'
  15. 'doublet_call'
  16. 'celltype.mapped_mnn'
  17. 'celltype.score_mnn'
  18. 'celltype_extended.mapped_mnn'
  19. 'celltype_extended.score_mnn'
  20. 'stage.mapped_mnn'
  21. 'cellstage.score_mnn'
  22. 'closest.cell_mnn'
  23. 'genotype'
  24. 'RNA_snn_res.1'
  25. 'seurat_clusters'
  26. 'across_genotype_integration_snn_res.1'
  27. 'extended_celltype'
  28. 'extended_celltype_mappingscore'
  29. 'original_celltype'
  30. 'original_celltype_mappingscore'
  31. 'stage_mappingscore'
  32. 'anatomy'
  33. 'anatomy_mappingscore'
  34. 'integrated_clusters'
  35. 'integrated_clusters_character'
  36. 'integrated_clusters_updated'
  37. 'integrated_clusters_updated_ordered'
In [85]:
DimPlot(object, 
        group.by = "integrated_clusters_updated_ordered") + theme(legend.position='false') + coord_fixed() + scale_color_manual(values=integrated_colors_updated)

Remove cells with high doublet scores¶

In [86]:
object$seurat_clusters_updated <- recode(object$seurat_clusters, 
  "12" = "1",
  "4" = "2",
  "10" = "3",
  "11" = "4",
  "6" = "5",
  "2" = "6",
  "5" = "7",
  "16" = "8",
  "17" = "9",
  "0" = "10",
  "9" = "11",
  "3" = "12",
  "15" = "13",
  "1" = "14",
  "8" = "15",
  "7" = "16",
  "13" = "17",
  "14" = "high_doublet_score",
  "18" = "low_quality")

# Assuming 'object' is your Seurat object containing the updated cluster information

# Define the desired order of levels
cluster_order <- c("0", "1", "2", "3", "4", "5", "6", "7", "8", "9", "10", "11", "12", "13", "14", "15", "16", "17", "high_doublet_score", "low_quality")

# Convert 'seurat_clusters_updated' to a factor with the specified levels and order
object$seurat_clusters_updated <- factor(object$seurat_clusters_updated, levels = cluster_order)
In [87]:
object_ss <- object[,!object@meta.data$seurat_clusters_updated %in% c("high_doublet_score", "low_quality")]
In [88]:
options(repr.plot.width=12, repr.plot.height=10)
DimPlot(object_ss, 
        group.by = "integrated_clusters_updated_ordered",
        split.by = "genotype") + coord_fixed() + scale_color_manual(values=integrated_colors_updated)
In [233]:
DefaultAssay(object_ss) <- "RNA"
options(repr.plot.width=12, repr.plot.height=10)
FeaturePlot(object_ss, 
        feature = "Hhex",
        order = TRUE,
        split.by = "genotype") + coord_fixed()
In [90]:
colnames(object_ss@meta.data)
  1. 'orig.ident'
  2. 'nCount_RNA'
  3. 'nFeature_RNA'
  4. 'barcode'
  5. 'sample'
  6. 'cell'
  7. 'mitochondrial_percent_RNA'
  8. 'ribosomal_percent_RNA'
  9. 'idx'
  10. 'stage'
  11. 'tdTom'
  12. 'tdTom_corr'
  13. 'pass_rnaQC'
  14. 'doublet_score'
  15. 'doublet_call'
  16. 'celltype.mapped_mnn'
  17. 'celltype.score_mnn'
  18. 'celltype_extended.mapped_mnn'
  19. 'celltype_extended.score_mnn'
  20. 'stage.mapped_mnn'
  21. 'cellstage.score_mnn'
  22. 'closest.cell_mnn'
  23. 'genotype'
  24. 'RNA_snn_res.1'
  25. 'seurat_clusters'
  26. 'across_genotype_integration_snn_res.1'
  27. 'extended_celltype'
  28. 'extended_celltype_mappingscore'
  29. 'original_celltype'
  30. 'original_celltype_mappingscore'
  31. 'stage_mappingscore'
  32. 'anatomy'
  33. 'anatomy_mappingscore'
  34. 'integrated_clusters'
  35. 'integrated_clusters_character'
  36. 'integrated_clusters_updated'
  37. 'integrated_clusters_updated_ordered'
  38. 'seurat_clusters_updated'
In [91]:
options(repr.plot.width=12, repr.plot.height=10)
DimPlot(object_ss, 
        group.by = "extended_celltype",
        split.by = "genotype") + coord_fixed() + scale_color_manual(values=extended.colors)

Generate Seperate WT and MT Seurat Object¶

In [79]:
table(object_ss@meta.data$genotype)
  mutant wildtype 
    4322     5648 
In [80]:
object_ss <- object_ss[,!object_ss@meta.data$integrated_clusters_updated_ordered %in% c("Low Quality")]
In [81]:
table(object_ss@meta.data$integrated_clusters_updated_ordered)
                     ExE Ectoderm                      ExE Endoderm 
                             1190                              2135 
                Visceral Endoderm                 Parietal Endoderm 
                             1420                                65 
                         Epiblast                  Primitive Streak 
                             1583                               797 
        Anterior Primitive Streak               Nascent mesoderm #1 
                              157                               721 
              Nascent mesoderm #2 Haematoendothelial progenitors #1 
                              491                                88 
Haematoendothelial progenitors #2  Posterior Lateral Plate Mesoderm 
                              169                               285 
                       Mesenchyme        Posterior Primitive Streak 
                               48                               173 
               Allantois Mesoderm             Intermediate Mesoderm 
                              523                                71 
                      Low Quality 
                                0 
In [82]:
object_ss_wt <- subset(object_ss, genotype %in% c("wildtype"))
In [83]:
object_ss_mt <- subset(object_ss, genotype %in% c("mutant"))
In [84]:
object_ss_wt
An object of class Seurat 
29998 features across 5614 samples within 2 assays 
Active assay: RNA (27998 features, 0 variable features)
 1 other assay present: across_genotype_integration
 2 dimensional reductions calculated: pca, umap
In [85]:
table(object_ss_wt@meta.data$integrated_clusters_updated_ordered)
                     ExE Ectoderm                      ExE Endoderm 
                              535                              1173 
                Visceral Endoderm                 Parietal Endoderm 
                              722                                41 
                         Epiblast                  Primitive Streak 
                             1569                               719 
        Anterior Primitive Streak               Nascent mesoderm #1 
                              156                               360 
              Nascent mesoderm #2 Haematoendothelial progenitors #1 
                              198                                28 
Haematoendothelial progenitors #2  Posterior Lateral Plate Mesoderm 
                                0                                 0 
                       Mesenchyme        Posterior Primitive Streak 
                                0                                96 
               Allantois Mesoderm             Intermediate Mesoderm 
                               17                                 0 
                      Low Quality 
                                0 
In [86]:
object_ss_mt
An object of class Seurat 
29998 features across 4302 samples within 2 assays 
Active assay: RNA (27998 features, 0 variable features)
 1 other assay present: across_genotype_integration
 2 dimensional reductions calculated: pca, umap
In [87]:
table(object_ss_mt@meta.data$integrated_clusters_updated_ordered)
                     ExE Ectoderm                      ExE Endoderm 
                              655                               962 
                Visceral Endoderm                 Parietal Endoderm 
                              698                                24 
                         Epiblast                  Primitive Streak 
                               14                                78 
        Anterior Primitive Streak               Nascent mesoderm #1 
                                1                               361 
              Nascent mesoderm #2 Haematoendothelial progenitors #1 
                              293                                60 
Haematoendothelial progenitors #2  Posterior Lateral Plate Mesoderm 
                              169                               285 
                       Mesenchyme        Posterior Primitive Streak 
                               48                                77 
               Allantois Mesoderm             Intermediate Mesoderm 
                              506                                71 
                      Low Quality 
                                0 

CellChat Pipleine¶

In [88]:
# Load required library
library(CellChat)

# Increase memory size limit (optional for large objects)
options(future.globals.maxSize = 50000 * 1024^2)

# --------------------------- #
# USER-DEFINED INPUTS
# --------------------------- #

# Define your Seurat object and set the cluster identity
input_object <- object_ss_wt
object_name <- "object_ss_wt"  # for file naming
Idents(input_object) <- "integrated_clusters_updated_ordered"

# Toggle downsampling (set to TRUE to enable)
enable_downsampling <- TRUE
downsample_n <- 4000
set.seed(123)

if (enable_downsampling) {
  cells.use <- sample(colnames(input_object), downsample_n)
  input_object <- subset(input_object, cells = cells.use)
}

# --------------------------- #
# PART 1: Create CellChat Object
# --------------------------- #

cellchat <- createCellChat(object = input_object, group.by = "ident", assay = "RNA")

# Load and assign database
CellChatDB <- CellChatDB.mouse
CellChatDB.use <- CellChatDB
cellchat@DB <- CellChatDB.use

# Show database categories (optional)
showDatabaseCategory(CellChatDB)

# Preprocessing
cellchat <- subsetData(cellchat)
# future::plan("multisession", workers = 4) # Uncomment for parallel processing
cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat)

# --------------------------- #
# PART 2: Inference
# --------------------------- #

cellchat <- computeCommunProb(cellchat, type = "triMean")
cellchat <- filterCommunication(cellchat, min.cells = 10)
cellchat <- computeCommunProbPathway(cellchat)
cellchat <- aggregateNet(cellchat)

# --------------------------- #
# PART 3: Visualizations
# --------------------------- #

# Define pathways to visualize
pathways.show <- c("NODAL", "WNT", "FGF", "BMP") 
netAnalysis_contribution(cellchat, signaling = pathways.show)

# --------------------------- #
# PART 4: Optional - Systems Analysis (Time-consuming)
# --------------------------- #

# Uncomment below if you wish to do functional/structural manifold analyses

# cellchat <- netAnalysis_computeCentrality(cellchat, slot.name = "netP")

# selectK(cellchat, pattern = "outgoing")
# cellchat <- identifyCommunicationPatterns(cellchat, pattern = "outgoing", k = 6)

# selectK(cellchat, pattern = "incoming")
# cellchat <- identifyCommunicationPatterns(cellchat, pattern = "incoming", k = 3)

# cellchat <- computeNetSimilarity(cellchat, type = "functional")
# cellchat <- netEmbedding(cellchat, type = "functional")
# cellchat <- netClustering(cellchat, type = "functional")

# cellchat <- computeNetSimilarity(cellchat, type = "structural")
# cellchat <- netEmbedding(cellchat, type = "structural")
# cellchat <- netClustering(cellchat, type = "structural")

# --------------------------- #
# PART 5: Save Output
# --------------------------- #

# Save with dynamic filename
output_dir <- "projects/10_Axin_1_2/github_submission/outputs/7/"
output_file <- paste0(output_dir, "cellchat_", object_name, ".rds")
saveRDS(cellchat, file = output_file)
[1] "Create a CellChat object from a Seurat object"
The `meta.data` slot in the Seurat object is used as cell meta information 
Set cell identities for the new CellChat object 
The cell groups used for CellChat analysis are  ExE Ectoderm ExE Endoderm Visceral Endoderm Parietal Endoderm Epiblast Primitive Streak Anterior Primitive Streak Nascent mesoderm #1 Nascent mesoderm #2 Haematoendothelial progenitors #1 Posterior Primitive Streak Allantois Mesoderm 
Issue identified!! Please check the official Gene Symbol of the following genes:  
 H2-BI H2-Ea-ps 
triMean is used for calculating the average gene expression per cell group. 
[1] ">>> Run CellChat on sc/snRNA-seq data <<< [2025-04-17 12:18:00]"
[1] ">>> CellChat inference is done. Parameter values are stored in `object@options$parameter` <<< [2025-04-17 12:20:18]"
In [89]:
# Load required library
library(CellChat)

# Increase memory size limit (optional for large objects)
options(future.globals.maxSize = 50000 * 1024^2)

# --------------------------- #
# USER-DEFINED INPUTS
# --------------------------- #

# Define your Seurat object and set the cluster identity
input_object <- object_ss_mt
object_name <- "object_ss_mt"  # for file naming
Idents(input_object) <- "integrated_clusters_updated_ordered"

# Toggle downsampling (set to TRUE to enable)
enable_downsampling <- TRUE
downsample_n <- 4000
set.seed(123)

if (enable_downsampling) {
  cells.use <- sample(colnames(input_object), downsample_n)
  input_object <- subset(input_object, cells = cells.use)
}

# --------------------------- #
# PART 1: Create CellChat Object
# --------------------------- #

cellchat <- createCellChat(object = input_object, group.by = "ident", assay = "RNA")

# Load and assign database
CellChatDB <- CellChatDB.mouse
CellChatDB.use <- CellChatDB
cellchat@DB <- CellChatDB.use

# Show database categories (optional)
showDatabaseCategory(CellChatDB)

# Preprocessing
cellchat <- subsetData(cellchat)
# future::plan("multisession", workers = 4) # Uncomment for parallel processing
cellchat <- identifyOverExpressedGenes(cellchat)
cellchat <- identifyOverExpressedInteractions(cellchat)

# --------------------------- #
# PART 2: Inference
# --------------------------- #

cellchat <- computeCommunProb(cellchat, type = "triMean")
cellchat <- filterCommunication(cellchat, min.cells = 10)
cellchat <- computeCommunProbPathway(cellchat)
cellchat <- aggregateNet(cellchat)

# --------------------------- #
# PART 3: Visualizations
# --------------------------- #

# Define pathways to visualize
pathways.show <- c("NODAL", "WNT", "FGF", "BMP") 
netAnalysis_contribution(cellchat, signaling = pathways.show)

# --------------------------- #
# PART 4: Optional - Systems Analysis (Time-consuming)
# --------------------------- #

# Uncomment below if you wish to do functional/structural manifold analyses

# cellchat <- netAnalysis_computeCentrality(cellchat, slot.name = "netP")

# selectK(cellchat, pattern = "outgoing")
# cellchat <- identifyCommunicationPatterns(cellchat, pattern = "outgoing", k = 6)

# selectK(cellchat, pattern = "incoming")
# cellchat <- identifyCommunicationPatterns(cellchat, pattern = "incoming", k = 3)

# cellchat <- computeNetSimilarity(cellchat, type = "functional")
# cellchat <- netEmbedding(cellchat, type = "functional")
# cellchat <- netClustering(cellchat, type = "functional")

# cellchat <- computeNetSimilarity(cellchat, type = "structural")
# cellchat <- netEmbedding(cellchat, type = "structural")
# cellchat <- netClustering(cellchat, type = "structural")

# --------------------------- #
# PART 5: Save Output
# --------------------------- #

# Save with dynamic filename
output_dir <- "projects/10_Axin_1_2/github_submission/outputs/7/"
output_file <- paste0(output_dir, "cellchat_", object_name, ".rds")
saveRDS(cellchat, file = output_file)
[1] "Create a CellChat object from a Seurat object"
The `meta.data` slot in the Seurat object is used as cell meta information 
Set cell identities for the new CellChat object 
The cell groups used for CellChat analysis are  ExE Ectoderm ExE Endoderm Visceral Endoderm Parietal Endoderm Epiblast Primitive Streak Anterior Primitive Streak Nascent mesoderm #1 Nascent mesoderm #2 Haematoendothelial progenitors #1 Haematoendothelial progenitors #2 Posterior Lateral Plate Mesoderm Mesenchyme Posterior Primitive Streak Allantois Mesoderm Intermediate Mesoderm 
Issue identified!! Please check the official Gene Symbol of the following genes:  
 H2-BI H2-Ea-ps 
triMean is used for calculating the average gene expression per cell group. 
[1] ">>> Run CellChat on sc/snRNA-seq data <<< [2025-04-17 12:20:44]"
[1] ">>> CellChat inference is done. Parameter values are stored in `object@options$parameter` <<< [2025-04-17 12:23:59]"
The cell-cell communication related with the following cell groups are excluded due to the few number of cells:  Anterior Primitive Streak 

CellChat Pipeline with Combined Datasets¶

In [12]:
## Step 1 ## Load in the Individual CellChat Objects

# set paths to data locations for the query samples
io$cellchat_wt = file.path(io$main, "projects/10_Axin_1_2/github_submission/outputs/7/cellchat_object_ss_wt.rds")
io$cellchat_mt = file.path(io$main, "projects/10_Axin_1_2/github_submission/outputs/7/cellchat_object_ss_mt.rds")

# load in the query data 
cellchat_wt <- readRDS(io$cellchat_wt)
cellchat_mt <- readRDS(io$cellchat_mt)
In [13]:
wt_celltypes <- unique(cellchat_wt@meta$ident)
mt_celltypes <- unique(cellchat_mt@meta$ident)

wt_mt_diff <- setdiff(wt_celltypes, mt_celltypes)
mt_wt_diff <- setdiff(mt_celltypes, wt_celltypes)

wt_mt_diff
mt_wt_diff
  1. 'Intermediate Mesoderm'
  2. 'Haematoendothelial progenitors #2'
  3. 'Posterior Lateral Plate Mesoderm'
  4. 'Mesenchyme'
In [14]:
# -----------------------------------------------
# Lift CellChat object to harmonize cell labels
# -----------------------------------------------

# Use the cell group labels from the more complete dataset (e.g., mt)
# This ensures both CellChat objects share a common set of cell group labels
group.new <- levels(cellchat_mt@idents)

# Apply the label structure to the WT object
cellchat_wt <- liftCellChat(cellchat_wt, group.new)

# Note: We do not apply liftCellChat to cellchat_mt, since it already contains all relevant cell groups

# -----------------------------------------------
# Merge CellChat objects for comparative analysis
# -----------------------------------------------

# Combine the lifted WT and original MT CellChat objects
object.list <- list(wt = cellchat_wt, mt = cellchat_mt)

# Compute centrality scores for each dataset
object.list$wt <- netAnalysis_computeCentrality(object.list$wt, slot.name = "netP")
object.list$mt <- netAnalysis_computeCentrality(object.list$mt, slot.name = "netP")

# Merge the objects, appending dataset names as prefixes to cell names
cellchat <- mergeCellChat(object.list, add.names = names(object.list), cell.prefix = TRUE)

# Note:
# - Warnings about barcode renaming are expected.
# - CellChat will reassign colnames in the expression matrix to match metadata.
# - This step merges data slots: 'meta', 'idents', 'net', 'netP', 'DB', etc.
The CellChat object will be lifted up using the cell labels ExE Ectoderm, ExE Endoderm, Visceral Endoderm, Parietal Endoderm, Epiblast, Primitive Streak, Anterior Primitive Streak, Nascent mesoderm #1, Nascent mesoderm #2, Haematoendothelial progenitors #1, Haematoendothelial progenitors #2, Posterior Lateral Plate Mesoderm, Mesenchyme, Posterior Primitive Streak, Allantois Mesoderm, Intermediate Mesoderm

Update slots object@net, object@netP, object@idents in a single dataset... 
Warning message in mergeCellChat(object.list, add.names = names(object.list), cell.prefix = TRUE):
“Prefix cell names!”
The cell barcodes in merged 'meta' is  LibB_WT_E6_1_IGO_14590_2#CTCCACATCCATCTCG-1 LibB_WT_E6_1_IGO_14590_2#CTCTCGAAGACGGTCA-1 LibB_WT_E6_1_IGO_14590_2#CGGCAGTAGTCCCTAA-1 LibB_WT_E6_1_IGO_14590_2#ACTATTCTCCGCAAAT-1 LibB_WT_E6_1_IGO_14590_2#TATACCTCAACAACAA-1 LibB_WT_E6_1_IGO_14590_2#GAGTTTGCAAGATGGC-1 
Warning message in mergeCellChat(object.list, add.names = names(object.list), cell.prefix = TRUE):
“The cell barcodes in merged 'meta' is different from those in the used data matrix.
              We now simply assign the colnames in the data matrix to the rownames of merged 'mata'!”
Merge the following slots: 'data.signaling','images','net', 'netP','meta', 'idents', 'var.features' , 'DB', and 'LR'.

In [15]:
cellchat
An object of class CellChat created from a merged object with multiple datasets 
 1028 signaling genes.
 8000 cells. 
CellChat analysis of single cell RNA-seq data! 

Identify Altered Interactions and Cell Populations¶

Compare the total number of interactions and interaction strength¶

In [16]:
options(repr.plot.width = 15, repr.plot.height = 6)

# Generate the first plot (count-based)
gg1 <- compareInteractions(cellchat, show.legend = FALSE, group = c(1, 2))
gg1$layers[[2]]$aes_params$size <- 5  # Increase number labels above bars
gg1 <- gg1 + theme(
  axis.text = element_text(size = 20),
  axis.title = element_text(size = 20),
  plot.title = element_text(size = 20, face = "bold"),
  legend.text = element_text(size = 20),
  legend.title = element_text(size = 20)
)

# Generate the second plot (weight-based)
gg2 <- compareInteractions(cellchat, show.legend = FALSE, group = c(1, 2), measure = "weight")
gg2$layers[[2]]$aes_params$size <- 5  # Increase number labels above bars
gg2 <- gg2 + theme(
  axis.text = element_text(size = 20),
  axis.title = element_text(size = 20),
  plot.title = element_text(size = 20, face = "bold"),
  legend.text = element_text(size = 20),
  legend.title = element_text(size = 20)
)

# Combine and display the two plots
gg1 + gg2
In [55]:
gg1 <- netVisual_heatmap(cellchat)
#> Do heatmap based on a merged object
gg2 <- netVisual_heatmap(cellchat, measure = "weight")
Do heatmap based on a merged object 


Do heatmap based on a merged object 


In [ ]:
help(netVisual_heatmap)
No documentation for ‘netVisual_heatmap’ in specified packages and libraries: you could try ‘??netVisual_heatmap’
In [57]:
options(repr.plot.width = 8, repr.plot.height = 8)
gg1
In [58]:
options(repr.plot.width = 8, repr.plot.height = 8)
gg2
In [17]:
options(repr.plot.width = 15, repr.plot.height = 6)

# Compare the major sources and targets in a 2D space
# (A) Identify cell populations with significant changes in sending or receiving signals
num.link <- sapply(object.list, function(x) {
  rowSums(x@net$count) + colSums(x@net$count) - diag(x@net$count)
})
weight.MinMax <- c(min(num.link), max(num.link)) # control the dot size in the different datasets

gg <- list()
for (i in 1:length(object.list)) {
  gg[[i]] <- netAnalysis_signalingRole_scatter(object.list[[i]],
                                               title = names(object.list)[i],
                                               weight.MinMax = weight.MinMax) +
    scale_color_manual(values = integrated_colors_updated) +
    theme(
      axis.text = element_text(size = 16),
      axis.title = element_text(size = 16),
      plot.title = element_text(size = 18, face = "bold"),
      legend.text = element_text(size = 14),
      legend.title = element_text(size = 16)
    )
}

patchwork::wrap_plots(plots = gg)
Signaling role analysis on the aggregated cell-cell communication network from all signaling pathways

Warning message:
“The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
of ggplot2 3.3.4.
ℹ The deprecated feature was likely used in the CellChat package.
  Please report the issue to the authors.”
Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
Signaling role analysis on the aggregated cell-cell communication network from all signaling pathways

Scale for colour is already present.
Adding another scale for colour, which will replace the existing scale.
In [45]:
options(repr.plot.width = 25, repr.plot.height = 7)

# (B) Identify the signaling changes of specific cell populations
gg1 <- netAnalysis_signalingChanges_scatter(cellchat, idents.use = "Nascent mesoderm #1", signaling.exclude = c("MK", "FN1"))
#> Visualizing differential outgoing and incoming signaling changes from NL to LS
#> The following `from` values were not present in `x`: 0
#> The following `from` values were not present in `x`: 0, -1
gg2 <- netAnalysis_signalingChanges_scatter(cellchat, idents.use = "Nascent mesoderm #2", signaling.exclude = c("MK", "FN1"))
#> Visualizing differential outgoing and incoming signaling changes from NL to LS
#> The following `from` values were not present in `x`: 0, 2
#> The following `from` values were not present in `x`: 0, -1
gg3 <- netAnalysis_signalingChanges_scatter(cellchat, idents.use = "Primitive Streak", signaling.exclude = c("MK", "FN1"))
#> Visualizing differential outgoing and incoming signaling changes from NL to LS
#> The following `from` values were not present in `x`: 0, 2
#> The following `from` values were not present in `x`: 0, -1
gg4 <- netAnalysis_signalingChanges_scatter(cellchat, idents.use = "Epiblast", signaling.exclude = c("MK", "FN1"))
#> Visualizing differential outgoing and incoming signaling changes from NL to LS
#> The following `from` values were not present in `x`: 0, 2
#> The following `from` values were not present in `x`: 0, -1
patchwork::wrap_plots(plots = list(gg1,gg2))
patchwork::wrap_plots(plots = list(gg3,gg4))
Visualizing differential outgoing and incoming signaling changes from wt to mt

Visualizing differential outgoing and incoming signaling changes from wt to mt

Visualizing differential outgoing and incoming signaling changes from wt to mt

Visualizing differential outgoing and incoming signaling changes from wt to mt

Warning message:
“ggrepel: 25 unlabeled data points (too many overlaps). Consider increasing max.overlaps”
Warning message:
“ggrepel: 27 unlabeled data points (too many overlaps). Consider increasing max.overlaps”
Warning message:
“ggrepel: 30 unlabeled data points (too many overlaps). Consider increasing max.overlaps”
Warning message:
“ggrepel: 26 unlabeled data points (too many overlaps). Consider increasing max.overlaps”

Part 2: Identify the altered signaling with network architecture and interaction strength¶

In [46]:
##cellchat <- computeNetSimilarityPairwise(cellchat, type = "functional")
#> Compute signaling network similarity for datasets 1 2
##cellchat <- netEmbedding(cellchat, type = "functional")
#> Manifold learning of the signaling networks for datasets 1 2
##cellchat <- netClustering(cellchat, type = "functional")
#> Classification learning of the signaling networks for datasets 1 2
# Visualization in 2D-space
##netVisual_embeddingPairwise(cellchat, type = "functional", label.size = 3.5)
#> 2D visualization of signaling networks from datasets 1 2
In [48]:
options(repr.plot.width = 15, repr.plot.height = 10)

# (A) Compare the overall information flow of each signaling pathway or ligand-receptor pair
gg1 <- rankNet(cellchat, mode = "comparison", measure = "weight", sources.use = NULL, targets.use = NULL, stacked = T, do.stat = TRUE)
gg2 <- rankNet(cellchat, mode = "comparison", measure = "weight", sources.use = NULL, targets.use = NULL, stacked = F, do.stat = TRUE)

gg1 + gg2
In [59]:
options(repr.plot.width = 15, repr.plot.height = 15)

i = 1
# combining all the identified signaling pathways from different datasets 
pathway.union <- union(object.list[[i]]@netP$pathways, object.list[[i+1]]@netP$pathways)
ht1 = netAnalysis_signalingRole_heatmap(object.list[[i]], pattern = "incoming", signaling = pathway.union, title = names(object.list)[i], width = 15, height = 15)
ht2 = netAnalysis_signalingRole_heatmap(object.list[[i+1]], pattern = "incoming", signaling = pathway.union, title = names(object.list)[i+1], width = 15, height = 15)
In [60]:
p1 <- draw(ht1, heatmap_legend_side = "right")
p2 <- draw(ht2, heatmap_legend_side = "right")
p2

Part IV: Visually Compare Cell Communication using Different Plots¶

In [61]:
pathways.show <- c("BMP") 
weight.max <- getMaxWeight(object.list, slot.name = c("netP"), attribute = pathways.show) # control the edge weights across different datasets
par(mfrow = c(1,2), xpd=TRUE)
for (i in 1:length(object.list)) {
  netVisual_aggregate(object.list[[i]], signaling = pathways.show, layout = "circle", edge.weight.max = weight.max[1], edge.width.max = 10, signaling.name = paste(pathways.show, names(object.list)[i]))
}
In [62]:
options(repr.plot.width=8, repr.plot.height=8)
pathways.show <- c("BMP") 
par(mfrow = c(1,2), xpd=TRUE)
ht <- list()
for (i in 1:length(object.list)) {
  ht[[i]] <- netVisual_heatmap(object.list[[i]], signaling = pathways.show, color.heatmap = "Reds",title.name = paste(pathways.show, "signaling ",names(object.list)[i]))
}
#> Do heatmap based on a single object 
#> 
#> Do heatmap based on a single object
ComplexHeatmap::draw(ht[[1]], ht_gap = unit(0.5, "cm"))
Do heatmap based on a single object 


Do heatmap based on a single object 


In [63]:
ComplexHeatmap::draw(ht[[2]], ht_gap = unit(0.5, "cm"))

Part V: Compare signaling gene expression distribtuion between different datasets¶

In [64]:
str(cellchat@meta$datasets)
 Factor w/ 2 levels "wt","mt": 1 1 1 1 1 1 1 1 1 1 ...
In [66]:
#cellchat@meta$datasets = factor(cellchat@meta$datasets, levels = c("NL", "LS")) # set factor level
options(repr.plot.width=10, repr.plot.height=12)
plotGeneExpression(cellchat, signaling = "BMP", split.by = "datasets", colors.ggplot = T, type = "violin")
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
In [67]:
#cellchat@meta$datasets = factor(cellchat@meta$datasets, levels = c("NL", "LS")) # set factor level
options(repr.plot.width=10, repr.plot.height=12)
plotGeneExpression(cellchat, signaling = "ncWNT", split.by = "datasets", colors.ggplot = T, type = "violin")
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
In [68]:
#cellchat@meta$datasets = factor(cellchat@meta$datasets, levels = c("NL", "LS")) # set factor level
options(repr.plot.width=10, repr.plot.height=12)
plotGeneExpression(cellchat, signaling = "WNT", split.by = "datasets", colors.ggplot = T, type = "violin")
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
In [150]:
#cellchat@meta$datasets = factor(cellchat@meta$datasets, levels = c("NL", "LS")) # set factor level
options(repr.plot.width=10, repr.plot.height=12)
plotGeneExpression(cellchat, signaling = "NODAL", split.by = "datasets", colors.ggplot = T, type = "violin")
There is no significant communication of NODAL

Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.
Scale for y is already present.
Adding another scale for y, which will replace the existing scale.

Load in the Axin Datasets Integrated with the Original Atlas¶

In [116]:
# set paths to data locations for the query samples
io$integrated_object = file.path(io$main, "projects/10_Axin_1_2/outputs/axin_object_integrated_e6_8_Original_Atlas.rds")

# load in the query data 
integrated_object <- readRDS(io$integrated_object)
In [125]:
io$integrated_object_clusters = file.path(io$main, "projects/10_Axin_1_2/outputs/integrated_object_clusters.rds")
integrated_object_clusters <- readRDS(io$integrated_object_clusters)
In [126]:
integrated_object <- AddMetaData(integrated_object, metadata = integrated_object_clusters)
In [127]:
colnames(integrated_object@meta.data)
  1. 'orig.ident'
  2. 'nCount_originalexp'
  3. 'nFeature_originalexp'
  4. 'cell'
  5. 'sample'
  6. 'embryo_version'
  7. 'stage'
  8. 'somite_count'
  9. 'anatomy'
  10. 'S_score'
  11. 'G2M_score'
  12. 'phase'
  13. 'louvain'
  14. 'leiden'
  15. 'celltype_PijuanSala2019'
  16. 'celltype_extended_atlas'
  17. 'sizeFactor'
  18. 'stage_updated'
  19. 'within_timepoint_integration_snn_res.2'
  20. 'seurat_clusters'
  21. 'originalexp_snn_res.2'
  22. 'within_atlas_version_integration_snn_res.2'
  23. 'nCount_within_timepoint_integration'
  24. 'nFeature_within_timepoint_integration'
  25. 'nCount_within_atlas_version_integration'
  26. 'nFeature_within_atlas_version_integration'
  27. 'nCount_RNA'
  28. 'nFeature_RNA'
  29. 'barcode'
  30. 'mitochondrial_percent_RNA'
  31. 'ribosomal_percent_RNA'
  32. 'idx'
  33. 'tdTom'
  34. 'tdTom_corr'
  35. 'pass_rnaQC'
  36. 'doublet_score'
  37. 'doublet_call'
  38. 'celltype.mapped_mnn'
  39. 'celltype.score_mnn'
  40. 'celltype_extended.mapped_mnn'
  41. 'celltype_extended.score_mnn'
  42. 'stage.mapped_mnn'
  43. 'cellstage.score_mnn'
  44. 'closest.cell_mnn'
  45. 'genotype'
  46. 'RNA_snn_res.1'
  47. 'across_genotype_integration_snn_res.1'
  48. 'integrated_clusters'
In [ ]:
table(object@meta.data$orig.ident)
Error in table(object@meta.data$orig.ident): object 'object' not found
Traceback:

1. table(object@meta.data$orig.ident)
In [140]:
object <- integrated_object

object <- object[,object@meta.data$orig.ident %in% "cell"]
In [141]:
object@meta.data$integrated_clusters_character <- as.character(object@meta.data$integrated_clusters)
In [142]:
object$integrated_clusters_updated <- recode(object$integrated_clusters_character, 
  "0" = "Epiblast",
  "1" = "ExE Ectoderm",
  "2" = "Epiblast",
  "3" = "Low Quality",
  "4" = "Nascent mesoderm #1",
  "5" = "Low Quality",
  "6" = "ExE Endoderm",
  "7" = "ExE Ectoderm",
  "8" = "Primitive Streak",
  "9" = "Allantois Mesoderm",
  "10" = "Intermediate Mesoderm",
  "11" = "Primitive Streak",
  "12" = "Low Quality",
  "13" = "ExE Endoderm",
  "14" = "Low Quality",
  "15" = "Primitive Streak",
  "16" = "Low Quality",
  "17" = "Low Quality",
  "18" = "Low Quality",
  "19" = "Posterior Lateral Plate Mesoderm",
  "20" = "Anterior Primitive Streak",
  "21" = "Mesenchyme",
  "22" = "Nascent mesoderm #2",
  "23" = "Visceral Endoderm",
  "24" = "Visceral Endoderm",
  "25" = "Low Quality",
  "26" = "Low Quality",
  "27" = "Epiblast",
  "28" = "Low Quality",
  "29" = "Haematoendothelial progenitors #2",
  "30" = "Haematoendothelial progenitors #1",
  "31" = "ExE Ectoderm",
  "32" = "Low Quality",
  "33" = "Low Quality",
  "34" = "Low Quality",
  "35" = "Posterior Primitive Streak",
  "36" = "Low Quality",
  "37" = "Low Quality",
  "38" = "Low Quality",
  "39" = "Low Quality",
  "40" = "Parietal Endoderm",
  "41" = "Low Quality",
  "42" = "Low Quality"
)
In [143]:
table(object@meta.data$integrated_clusters_updated)
               Allantois Mesoderm         Anterior Primitive Streak 
                             2219                              1611 
                         Epiblast                      ExE Ectoderm 
                             9337                              7860 
                     ExE Endoderm Haematoendothelial progenitors #1 
                             3535                              1220 
Haematoendothelial progenitors #2             Intermediate Mesoderm 
                             1191                              2624 
                      Low Quality                        Mesenchyme 
                            28863                              1707 
              Nascent mesoderm #1               Nascent mesoderm #2 
                             2651                              1236 
                Parietal Endoderm  Posterior Lateral Plate Mesoderm 
                              211                              1558 
       Posterior Primitive Streak                  Primitive Streak 
                              802                              7093 
                Visceral Endoderm 
                             1841 
In [144]:
# Assuming 'object' is your Seurat object containing the updated cluster information
# Define the desired order of levels
cluster_order <- c("ExE Ectoderm", "ExE Endoderm", "Visceral Endoderm", "Parietal Endoderm",
                   "Epiblast", "Primitive Streak", "Anterior Primitive Streak", "Nascent mesoderm #1", "Nascent mesoderm #2",
                   "Haematoendothelial progenitors #1", "Haematoendothelial progenitors #2", 
                   "Posterior Lateral Plate Mesoderm", "Mesenchyme", "Posterior Primitive Streak", "Allantois Mesoderm", "Intermediate Mesoderm", "Low Quality")
In [145]:
# Convert 'seurat_clusters_updated' to a factor with the specified levels and order
object$integrated_clusters_updated_ordered <- factor(object$integrated_clusters_updated, levels = cluster_order)
In [147]:
DimPlot(object, reduction = "umap",
        group.by = "integrated_clusters_updated",
        split.by = "stage") + theme(legend.position='false') + coord_fixed() + scale_color_manual(values=integrated_colors_updated)
In [212]:
DimPlot(object, reduction = "umap",
        group.by = "integrated_clusters_updated") + theme(legend.position='false') + coord_fixed() + scale_color_manual(values=integrated_colors_updated)
In [236]:
DimPlot(object, reduction = "umap",
        group.by = "celltype_PijuanSala2019") + theme(legend.position='right') + coord_fixed() + scale_color_manual(values=original.colors)
In [238]:
options(repr.plot.width=24, repr.plot.height=10)

DimPlot(object, reduction = "umap",
        group.by = "celltype_extended_atlas") + theme(legend.position='right') + coord_fixed() + scale_color_manual(values=extended.colors)
In [152]:
DefaultAssay(object) <- "originalexp"
In [241]:
FeaturePlot(object, reduction = "umap",
        feature = "Lhx1",
        order = TRUE) + coord_fixed()
In [245]:
DotPlot(object, group.by = "celltype_PijuanSala2019",
        feature = c("Lhx1", "Sox17", "Sfrp1", "Gpc4", "Otx2", "Chrd", "Foxa2", "Cer1"))
In [ ]: